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This expansion is particularly useful when considering ionospheric propagation at low 
frequencies. The complex problem dealing with two media, viz., a homogeneous earth 
and a surrounding stratified atmosphere, leads to intractable expressions. However, as 
the influence of the earth may be accounted for by an approximate boundary condition at 
the earth's surface, the problem is then reduced to that of the outer medium only. The 
coefficients of the mode expansion for this simplified problem will be derived while taking 
into account the earth's curvature; however, the latter proves to be negligible under very gen- 
eral conditions. The expansion to be derived is wanted in particular when studying the 
influence of a gradual transition in the electron density with height at the lower edge of the 
ionosphere. 

1. Introduction 

Maxwell's equations for time-harmonic functions proportional to exp (—iut) can be put 
in the following form (expressed in Gaussian units) provided the dielectric constant e, the 
permeability \x and the conductivity a arc independent of time: 



curl e—i-iJih=0', curl h+i- e c 6=0, 
div ye c e)=0' } div (^)=0; 



(1) 



the complex parameter e c =e-\-i4ir2/u includes both dielectric and conductive properties. 
We take /x=l and assume a stratified atmosphere of spherical symmetry, all properties 
of which can be considered as functions of the distance r to the center of the earth. We ac- 
cordingly substitute € c =k 2 (r)jk\, k =o)/c being the wave number in a vacuum; k(r) can be 
interpreted as the local wave number at a height r— a (a = earth's radius) above the surface 
of the earth. In the absence of the earth's magnetic field e c (r) constitutes a scalar. The 
Maxwell equations (1) can then be satisfied by the field 



= , 2 , , curl curl < k(r)Hr >;A=r-curl < k(r)Hr >> 



(2) 



W 



provided the Hertzian vector H r has a radial direction throughout space (r being a vector of 
length r in this direction) ; moreover, the amplitude H has to satisfy the following scalar wave 
equation in the outer space r^>a: 

(A+k 2 eS )H=0, (3) 

in which the new parameter k 2 eS is defined as follows [l] 1 : 

ftW -* W -*(r)£{l^}. 

The solution (2) is of the so-called electric type which includes the case of a vertical electric 
dipole to be considered hereafter. Let the dipole be placed in a homogeneous slab extending 
in the region a<V<&' directly above the earth's surface, the upper boundary r=a' of which 
may be identified with the lower edge of the ionosphere in the case of ionospheric propagation. 

* Present address: Philips Research Laboratories, Eindhoven, The Netherlands. 
'Figures in brackets indicate the literature references at the end of this paper. 
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Assuming a wave number k inside this slab, the primary field of an infinitesimal vertical dipole 
can be derived from (2) by substituting for H the function: 

„ p /Z exp ftTP) , , . 

T and P here represent the positions of transmitter and receiver, TP their mutual distance, b 

the distance from T to the center of the earth, and // the moment of the transmitter. The 

-> 

correctness of (5) is shown by comparing the associated radial Hertzian vector H vr r with the 
conventional Hertzian vector for such a dipole. The latter Hertzian vector has a constant 
direction (parallel to the dipole axis) throughout space, its amplitude being given by 

II exp (ihTP) 
c TP 

It can be verified that both Hertzian vectors only differ by a gradient, namely that of the 
quantity ill exp(ik TP) / (bk c) ; however, such a gradient has no effect whatever on the primary 

field which can therefore also be derived from the radial Hertzian vector H pr r according to (5). 

2. Form of the Mode Expansion 

The homogeneous wave equation (3) does not hold at the transmitter itself. The latter 
involves a singularity of the amplitude H at the point T. In view of (5) this singularity is 
the same as that of the function Il/(cbTP). We introduce polar coordinates with the origin 
at the center of the earth and the axis 0=0 drawn through the transmitter, the position of 
which therefore is given by r=b, 0=0. The singularity in question can be accounted for by 
a right-hand side of (3) which only differs from zero at T. The complete wave equation in 
the outer space r^>a thus proves to be: 

(A+ki s )H=-±^8(b-r) 8 -f- (6) 

This equation can be verified by integrating both sides over an infinitesimal sphere around T, 
while applying Gauss's integral theorem to the left-hand side, and the following properties of 
the impulse function to the right-hand side: 

P°^-sin0d0= ( $ ° d{ff)dd=l; (^ 8(b-r)dr=l. 

JO Jo * Jft-e 

We next assume that the influence of the earth may be accounted for by the boundary 
condition 

-^(rH)=yrH, (r=a) (7) 

which is the equivalent of the corresponding condition ()H/c)z=yH for a plane boundary. The 
quantity — iy/{koe c (a) } can be interpreted as a surface admittance (ratio at grazing incidence 
of the electric field and the magnetic field, both at infinitesimal distances above and below 
the earth's surface). For e c (a) = l the parameter y is given, for a vertical dipole, by 

&! being the complex wave number inside the earth. The special case of a perfectly conducting 
earth corresponds to 7=0. The value of 7 can also be chosen such as to account for a strati- 
fication in the ground [2]. 
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The approximative boundary condition (7) has two advantages in the present problem, 
viz., (a) it reduces the original two-media problem to that of the outer medium only, (b) it 
guarantees the existence of a complete set of discrete orthogonal modes in the outer space 
a<><°°. These modes are characterized as solutions of the wave equation (3) that satisfy 
both the boundary condition (7) at r=a and a proper radiation condition at r= oo . 

The general form of the modes is obtained by determining solutions of (3) with the aid 
of a separation of variables. Owing to the axial symmetry of the present problem around 
the axis 6=0 (containing the transmitter) we only have to concern the polar coordinates 
r and 0. The expansion of the scalar H in terms of the modes in question then becomes of 
the type: 

H=Jbcif ni (r)P ni { cos (t-«) } ; (8) 

the height-gain functions f n (r) occurring here are solutions vanishing at infinity of the differ- 
ential equation 



£ (ftt + {tt(r)-2^}w„)-0; 



(9) 



finally the discrete eigenvalues Ui (constituting the order of a Legendre function) are to be 
determined from the boundary condition 

j r {rf n )=yrf n {r=a). (10) 

3. A First Evaluation of the Coefficients of the Mode Expansion 

The condition (10) guarantees the required boundary condition (7) for the complete solution 
(8). All Legendre functions in (8) become infinite along the whole axis 6=0, n x being not an 
integer; however, the proper choice of the coefficients c, x will reduce this singularity to the 
required singularity at the transmitter T only. In view of properties of Legendre functions 
the singularity of each term of (8) appears from the approximation 

9 R 2 n 

P n { cos (tt—6) } ~— sin (mr) log -=- sin (nw) log -— > (0->O or p->0) (11) 

TV Ji T LV 

p being the distance from the point (r, 0) to the axis 0=0. The approximation (11) can be 
proved by applying Tauber's theorem [3] 



Lim< pf(p) l=Lim-^ 



to the operational relation h(t)=f(p) for which [3, p. 243, eq (35)] 

r(p)r( P +i) 



h(t)=P n (2e- i -l)U(t); f(p)- 



T(p+n+l)T(p-ny 



while substituting t=—2 log sin (0/2). In its turn the approximation (11) involves the three- 
dimensional equation : 

CM-ifc) [Ur)P n {cos (r-0) }]=l sin M M ^L (12) 

This equation can be proved along the same lines as (6), but by integrating over an infinitesimal 
cylinder around 0=0 instead of an infinitesimal sphere around T; the proof is facilitated by 
replacing 5(0)/(r 2 0) by the equivalent quantity 8(p)/p. 
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A summation of (12) over all eigenvalues, after an initial multiplication by c n , yields, in 
view of (8) : 

(A+*k) H=l ~ S -^-± Cl sin MOAW- 

n I U 1= — oo 

The introduction of both negative and positive modes avoids the occurrence of a deviation 
of c from the general expression for c x valid for 1^0. In fact, it is well known from the special 
case of a perfectly conducting earth and ionosphere that the addition of the I and — I mode to 
a single one involves an exceptional value of c (Neumann factor). Obviously, the right-hand 
side of the latter relation has to be identical with that of (6), so as to have: 

1 CD J] 

—2 Jlji sin (im,)/ ff| (r) = -^ 8(b -r) . 
A further multiplication by r 2 yields: 

~ ± c t sin (Tn t )f ni (r)=~8(b--r). (13) 

7T Z=_co l CO 

The right-hand side of (13) can also be expanded in terms of the functions /"^ (r). In fact, 
the complete orthogonality of the latter over the interval a<jr<^ o° involves the following 
expansion for an arbitrary function (p(r) : 



f* 00 

<p(r)f ni (f)dr 
«>(»■)= S ^~ -fn t (r). 



An application to (p(r-) = b(b—r) results in: 

* =_ " fl,(r)dr 

J a 

Hence, by equating the coefficient of f n (r) in both sides of (13), we can evaluate c x . The 
substitution of the value thus found into (8) yields the final form of the mode expansion, viz., 

„ II ^ fn,(b)fn,(r) P.,{COB fr-g)} 

Jti = —ir -T- 2-? ~7To3 r^ v (14) 

cb i^Too I j>2f\-, (sin 7m f ) 

f ni (r)dr 

J a 

4. Another Representation of the Mode Expansion 

The coefficients of (14) only depend on the eigen functions/^ at the transmitter and receiver, 
apart from the integral in the denominator. Fortunately, the latter can be reduced as follows 
to a quantity referring to the earth's surface; hence, the behavior of j n elsewhere needs not 
to be known explicitly. We apply the operator 

Ijfr/O-l-r/.-l; 

to (9). The resulting expression can be represented in the form: 
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in which d/dr denotes a differentiation with respect to r at constant n. An integration of 
this equation over the interval a<jr<i », taking into account the vanishing at infinity of the 
expression between braces, yields a relation which can be put in the form: 



(2n+l)J"f n (r) tfr=- r 2 f 2 ^J dr (r/J L 



We bear in mind that (d/dr) (rf n )/(rf n ) at r=a does in general depend on n, though this 
quantity assumes the value 7 independent of n% at the eigenvalues n x of n [see (10)]. The 
integral in (14) can thus be evaluated; it leads to the following alternative form of the mode 
expansion for H: 



ca 2 b i^r m 



(2n,+l) 


/», (&) /«, (f)P ni { cos (x-0)} 


L I r /« J J 


/2, (a) sin (im,) 



(15) 



We now pass to the vertical component E r of the electric field, the most important field 
component for the vertical dipole under consideration. The field representation (2), com- 
bined with the wave equation (3), implies the formula: 



Er=i(v„+^m. 



The operator k 2 ea -\-b 2 /dr 2 only affects the height-gain factors /^(r) of (15). The resulting 
quantity reads, in view of (9), 



J n m 



We thus obtain from (15) the following further mode expansion: 

iirll a n t (rii + l) (2n t + 1) fn l (b)fn l (r) Pn,{cos (tt-0)} 



E r - 



kca 2 br l== _ o 



-1 dr 



d I 4- (rfj 



dn j 



rf„ 



f 2 ni (a) 



sin (xni) 



(16) 



r — a; n=nj 



5. Connection of the Mode Amplitude With the Atmospheric Reflection 

Coefficient 

The dependence of the mode coefficients in (16) on the height-gain functions f n (r) can be 
replaced by that on a properly defined reflection coefficient. In order to discuss this coefficient 
we assume a homogeneous layer a<Cr<Cjob' with wave number k=k above the earth's surface. 
In the case of ionospheric propagation this layer may be identified with the space between 
the earth and the ionosphere, its upper boundary r=a' then representing the lower edge of 
the ionosphere (at a height h = a'—a above the earth's surface). In the case of the troposphere 
the thickness h might finally approach zero; the wave number in the homogeneous slab then 
had better be taken equal to k' 5*k Q in order to account for the small difference of the refractive 
index at the earth's surface from its vacuum value of unity. We shall restrict ourselves to 
the ionospheric case; it then is sufficiently accurate to identify the wave number in the inter- 
space between the earth and the ionosphere with its vacuum value k . 

In the homogeneous slab the general solution of (9) is well known, viz., a combination 
of two spherical Hankel functions: 



(1) / tt \H 



(1) 



79 



(1) 

,(2)/ 



The corresponding complete wave functions h™ (k Q r)P n {cos (tt—6)} represent a rising and a 
descending wave, respectively. A special mode will in general contain both these waves. 
Hence the height-gain function has the form: 



Mr) =ajii? (k r) +p n h™ (fcor) (a<r<a') . 



(17) 



Near the earth's surface the logarithmic derivative in the denominator of (16) therefore can 
be represented as follows: 



J: (rfn) cc n | {rft™ (k r) } +0 n | {rhi« (for) } 



Kf. 



r{a n h^(k r)+fi n h^(k r)} 



(18) 



This quantity depends, amongst others, on the ratio a n /P n which can be connected with 
a reflection coefficient. In fact, the second term of (17) constitutes the reflected wave that 
corresponds to the rising wave represented by the first term; hence we can interpret the ratio 
of both terms at r—a' as an atmospheric reflection coefficient T(n). We then obtain: 



1{n) ~a n h^(k a')' 



(19) 



with the aid of which (18) can be transformed as follows at r=a: 



rfn 



i{rhlP(k r) } T=a +T(n) ^%?,\ • i{rh£> (fr) } r =a 



dr 



h n 2, (k a') dr x 



»{ 



W\ha)+T(n) |^|^| Ai 2) M \ 



(20) 



For a stratification degenerating to a homogeneous medium in the domain r>a' the coefficient 
T(n) proves to reduce to the Fresnel reflection coefficient for an elevation angle (compare next 
section) r=arc sin (n/k Q a) if, moreover, the earth's curvature is neglected. 

The rigorous expression (20) has to be differentiated with respect to n in order to obtain 
the denominator of (16). We shall derive a suitable approximation of this complicated quan- 
tity in the next section; it is then convenient to convert (20) into the alternative form: 



rfn 



{rh^(k r)} 



-. 1 + T(n) 



rh^ihr) 






{rh^ihr). 



d 



f r {rh^(hr)\ 



l + T(n) 



te l >(JU')fc 2) 0U) 
hn^ikoa^h^ikoa) 



(21) 



6. An Approximation of the Mode Expansion (16) 

The approximations for the function h„(z) have a different analytical form according to 
the position in the complex plane of the parameter z/n. We restrict ourselves to those values 
of z/n that correspond (for real z~£ eO to eigenvalues n t of n of type I, using Rydbeck's no- 
menclature [4]. These eigenvalues can be characterized by the unequalities : 

~Ren<.k a— (ha) 1 ' 3 , lm»>0>|»l. 

Such eigenvalues prove to play the dominating role in long-wave propagation. The approxi- 
mations valid for the corresponding range of z/n read [5] : 
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exp|(-/,.+ i) r n (l-« 2 ) 1/2 vi ex P 

A a» (2) ^__L. J» 1/4 M J ; A< 2 >( 2 )~i — 



-(^+i)J i (i-u 2 ) 1/2 



" 



I 



2\l/4 



(22) 



in which the square roots are assumed as having a positive real part. 

In view of the quantities occurring in (21) we are particularly interested in the two 
following further approximations which can be derived at once from (22): 



||]^^exp{2 ( . + i ) j; / \l-. 2) v 2 ^, 



d_ 
dz 



zh${z) 



1 , s 
^2z~ r 2(n 2 -z 2 ) 



2~2\l/2 



(n 2 -z 2 ) 



(23) 



(24) 



We next substitute n=k a sin r or n=k a' sin r'. The asymptotic approximations of 

(i) 
the complete wave functions h™ (k r) P re {cos (x— 0)} then show that r and r' can be in- 
terpreted as the (in general complex) angles between the ray trajectories associated with the 
mode in question, and the verticals at r=a and r=a' respectively. Moreover, we substitute 
for z k Q a or k a' . The corresponding transformations of (23) and (24) lead, amongst others, 
to the following approximations: 



h™(k a 

K 2) (k a) 



Ty^iexp 2ik < a/ cos r' — a sin r ( ~ r )\ \ 
xp —2ik () a J cos r— sin r(^—r J r \, 



~% ex 



'&{«*>«*■)} 



(i) 



rh^(hr) 



-*b( 



Ti cos r- 



tan 2 1 



> 



An application of the two latter relations yield : 



J WW); 



|;WW); 







^ K 2) (k a) 
A h^(k a) 



7 /. tan 2 . . 

-k l I COS r — — rrz - } 



COS T- 



— J i exp — 2ifc aJ cos t— sin r (o"~t Iff 



The final substitutions of all these approximations into (21) yields: 
d 



< 



dr 



iTfn) 



rf n 



j — (ik 



COS T- 



tan 2 T 
~~2a 



> 



1 _ T(n) e 2i *o" «»" r l+it^h/(2koacosT) 
1—i tan 2 T/(2& a cos r) 



(25) 



l + T(ri)e 2i W* 



provided we neglect the small difference between r and r' (this difTerence vanishes altogether 
in the flat-earth approximation). 
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In (16) we need the derivative of (25) with respect to n. According to the relation n 
=k Q a sin r we have 

j^ = 1_ j^ 

dn k a cos r dr 

A corresponding evaluation of the ^-derivative of (25) results in: 









1 f . . sin r / . , tan 2 r\ d ^ 

= —— -Kt sm r— y— — * — ( i cos t+-^ — )^— v 
a cos r ^ k a cos 6 r \ 2k a Jot J 



k a cos 3 r V" w ° ' ' 2fc a 

1 _ Tin) e 2ik o h cos - 1+ ^ tan V( 2 ^ CQS r ) 
1— i tan 2 r/(2/:oa cos r) 



(26) 



l + T(^)6 2 ^0^os T 



Moreover, we shall apply the well-known asymptotic approximation for Legendre func- 
tions having a large complex order n with positive imaginary part, viz., 



P n { COS fr-fl)} 

sin (t27t) 



/_2|_y* e 

\7T7i sin 0/ 



«*(»+ 1/2)0 



(27) 



A final substitution of both (25) and (27) into (16), while approximating rii+1 and r^+1/2 
by ni=k a sin r*, yields the following expansion in the simplest case of a transmitter and 
receiver on the ground (r=b=a): 



c \a sm 6/ 



sin 5/2 r; e ik Q ad sinr * 



{ % tanr / i sin 2 r \ d \ 

k a cos 3 r \ 2t a cos 3 r/dr J 



,3,^ 



' i _ zv n ) 6 2i*b» co Sr l+^sin 2 r/2fc acos 3 r 

_____ l—i sm 2 r/2k a cos 3 r 

l + T(n)e 2ik o hcOBT 



(28) 



7. Flat Earth Approximation of the Mode Expansion 

The following simplification results from (28) for a= oo } remembering the relation 6=d/a 
(d= distance from transmitter to receiver) : 



J? _9 k /211 ( 2iri \ n ^ S j n 5/2 r?6 ^sln r; 

^ T ^ c \d ) i~«r/ t d \C \-T{n)e 2i W^ \-\ ' 

|_\ T dJ\l + T(7l)6 2 ^ cos r j- J 



(29) 



This expression can also be derived straightforwardly by applying the flat earth condi- 
tions right from the beginning, while replacing the boundary condition (7) at the earth's 
surface by 

|L= 7 n (2=0). (30) 

The rigorous expansion for the amplitude of the vertically directed Hertzian vector II [to be 
substituted in (2) instead of Hr ] then arrived at reads: 



c 



i^.p fdh/dz\l 



A,(Wx,(s) 

A 2 ,(0) 



B5 U (M). 



(3D 



z=0;X=X, 
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In this series, replacing (15), f\(z) denotes a height-gain function solving the differential 
equation 

^+{/4(2)-x 2 }/a=o, 

and satisfying the radiation condition at infinity (z= »). The other boundary condition (30) 
involves the relation 

A'(0)=7/x(0) (0=0) (32) 

at the earth's surface [replacing (7)]; this condition can only be satisfied if X equals one of the 
eigenvalues X?; h represents the height of the transmitter above the earth. The expansion 
for the vertical component of the electric field [replacing (16)], to be derived from (31), 
becomes : 

,_ 2?r n^ M fx^f^jz) 

k c i=-oo \o_ J dh/dz \ Jx l (0) 



ra fdh/dz\i 

U>X\ h /Ji-osx-x, 



The wave functions in the homogeneous slab o^z^h above the earth here reduce to 
simple exponentials so that (17) is to be replaced by: 

f x (z) =a x cxp { - (\*-k 2 )^z } + ft exp { (\*-k 2 ) l "s } , (34) 

which enables an easy evaluation of the denominator in (29). The final series representing e z 
for a transmitter and receiver on the ground (h o =z=0) then proves to be, substituting more- 
over \i=k s'u\ T t , 

e z = —^ IlJ£ m ry— ^\f i-T(r)e^ C0 ^ \-] " //(J) (W Sin Tz) * (35) 

This rigorous series leads to (29) if the Hankel function is replaced by its asymptotic approxi- 
mation. A comparison of (29) and (28) shows under what circumstances the earth's curvature 
becomes negligible. 

An evaluation of the derivative in the denominator of (35) leads to the alternative form: 

2irikl TJ ^S sin 3 r z H^(k d sin r t ) 



I tan T ]i y^ifco* cos T ^ Q I Tflityi cos r \2 

The denominator of (36) can further be worked out by considering the equation determining 
the modes. The latter is arrived at as follows. A substitution of (34) while evaluating the 
boundary condition (32) leads to the relation 

gfeO=l(Xi-4i)i/i_ 7 £opX=X,- 

{p/a)+l 

On the other hand, the definition of the atmospheric reflection coefficient for the flat case 
involves the following analogy of (19) : 

a 
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An elimination of pfa from the two latter relations yields an equation which can be put in the 
following form [replacing X by k smr and (X 2 — &o) 1/2 by —ik cos r]: 

ik cos Tj+y ^ T ,^ e 2 iko hco8T l= i 
ik cos T Z —7 

The first quotient here proves to be nothing else but the approximative plane-wave reflection 
coefficient at the earth's surface, R{rj) say, that results from the boundary condition (32). 
Hence the mode equation reads : 

R(r l )T(r l )e 2ik o hcoaT i=l. (37) 

This simple equation is also obtained at once from the following resonance condition. We 
consider a plane rising wave in the homogeneous slab 0<£</i, and derive the other plane rising 
wave that results from it after a reflection against the lower edge of the stratification and a 
reflection thereafter against the earth surface. The relation (37) then expresses the condition 
that the latter rising wave should be identical to the original one. 

The exponential exp (2ikoh cos t t ) in (36) can now be replaced throughout by {R(ti) T(r z )} _1 . 
This yields the following expansion replacing (36) : 

2_«§ £ { ]i,Z£i ] l Sin3 \, Tr> , - m <M sin rO • (38) 

c i±r m {(1 — R 2 ) tan r+2R(2ik h suit— dlog T/or} T=Ti 

For completeness sake we also mention the corresponding asymptotic approximation based 
on large values of k d, viz., 

(Ti)K(2fc )* n ^ {l+flWl'flin^r.cosT, 1W8ta . 

* C dH '--- { l-g+2fi cos x (2tU- » l0 | n ? T ) } T _ r; ' 

8. Special Case of a Flat Stratification Reducing to a Homogeneous 

Atmosphere 

A further reduction of (38) is only possible by introducing a special model for the stratification 
in the space z^>h. Let us consider the example of a homogeneous atmosphere, sharply bounded 
at z=h. The other boundary z=0 of the interspace 0<2</i, also being sharp, the situation 
becomes almost identical to that of a wave guide limited by two infinite plane parallel plates 
with a separation h. However, we have to bear in mind the approximative boundary condition 
applied at 2=0 so that this wall behaves differently from the other wall &tz=h. Let us mark the 
refractive index of the atmosphere by %\\ the reflection coefficient T then reduces to an ex- 
pression of the Fresnel type, viz., 

j,, , n 2 cos r— (n 2 t — sin 2 T)K 

n 2 cos r + (n 2 — sin 2 r) 1/2 

This formula involves the relation 

dlogT 7 ^ 2tanr(l + r) 
dr (1 — T)(sin 2 r— n 2 cos 2 r)' 

with the aid of which (38) can further be transformed [6] into : 

(1 + i?) 2 sin 2 r cost 



C Z=-oo 



{1 + T *\ 

ikjh COST — -; ^ w . g = rr V 
(l — T) (sinV— n 2 cos 2 r) J 



H ( ?(MsinT,). (40) 
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Owing to the approximative boundary condition at 2—0 we here obtain a completely discrete 
The rigorous boundary condition at z=0, on the contrary, would yield also an integral 
contribution which can be ascribed to a continuous part of the mode spectrum. However, the 
latter vanishes in the case of a perfectly conducting wall at z=0. 

The result (40) is also rigorous for a wave guide limited by two perfectly conducting walls at 
2=0 and z = h;we then have T=R=1. In this case we derive from (40) the following series for 
the vertical field component at z=0 that is due to a dipole also situated at z=0: 

27T&Q 



II SsiiA-jff^ (MsinT,). 

Clb Z=-oo 



The corresponding asymptotic approximation becomes : 
moreover, we may substitute r^=arc cos (irl/koh). 
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